# Filtros gaussiano y bilateral.

import nibabel

from pylab import *

import scipy.spatial.distance
import scipy.stats
	
def gaussianFilter(img, sigma, z):	
    #3 sigma para cada lado, mas uno para que sea impar
    window_size = int(6. * sigma)
    
    if window_size % 2 == 0:
        window_size += 1
    
    mid = window_size / 2
    
    #print img
    #print mid, window_size
    
    #Creamos la variable correspondiente a la ventana, con tipo float, sin inicializar los valores
    window = empty((window_size, window_size, window_size), dtype = float)
    window[mid, mid, mid] = 1. 
    #print window[mid, mid, mid]
    
    for i in xrange(window_size):
        for j in xrange(window_size):
            for k in xrange(window_size):
                if(i == mid and j == mid and k == mid):
                    continue
                
                #value = max(abs(i - mid), abs(j - mid), abs(k - mid))
                value = scipy.spatial.distance.euclidean([i, j, k], [mid, mid, mid])
                window[i, j, k] = scipy.stats.norm.pdf(value, 0, sigma)

    #print window    
    
    l = empty((img.shape[0], img.shape[1], img.shape[2]), dtype = float)
	
    for i in xrange(img.shape[0]):    
        for j in xrange(img.shape[1]):
            for k in xrange(z, z + 1):
                l[i, j, k] = conv(img, window, i, j, k)
                
    #print l
    
    return l    

def bilateralFilter(img, sigma_1, sigma_2, z):	
    #3 sigma para cada lado, mas uno para que sea impar
    window_size = int(6*sigma_1)
    
    if window_size % 2 == 0:
        window_size += 1
    
    mid = window_size / 2
    
    #print window_size, mid
    
    #Creamos la variable correspondiente a la ventana, con tipo float, sin inicializar los valores
    window = empty((window_size, window_size, window_size), dtype = float)
    window[mid, mid, mid] = 1.    
    
    for i in xrange(window_size):
        for j in xrange(window_size):
            for k in xrange(window_size):
                if(i == mid and j == mid and k == mid):
                    continue
                
                value = max(abs(i - mid), abs(j - mid), abs(k - mid))
                window[i, j, k] = normpdf(value, 0, sigma_1)
                
    #print window
    
    l = zeros((img.shape[0], img.shape[1], img.shape[2]), dtype = float)    
    
    for i in xrange(mid, img.shape[0] - mid):
        for j in xrange(mid, img.shape[1] - mid):
            for k in xrange(z, z + 1):
                window_intensity_local = calculateIntensityMatrix(getSubmatrix(img, [i, j, k], [mid, mid, mid]), sigma_2)
                #print window_intensity_local                
                window_bilateral = multiply(window, window_intensity_local)
                l[i, j, k] = conv(img, window_bilateral, i, j, k)
    return l

# Calculo la matriz para toda la imagen
def calculateIntensityMatrix(img, sigma):
    window_size = img.shape
    mid = []
    
    for i in xrange(3):
        mid.append(window_size[i] / 2)
        
    #print mid
    
    # Ventana para la diferencia de colores
    window_intensity = empty((window_size[0], window_size[1], window_size[2]), dtype = float)    
   
    for i in xrange(window_size[0]):
        for j in xrange(window_size[1]):
            for k in xrange(window_size[2]):
                value = abs(img[i, j, k] - 1.0 * img[mid[0], mid[1], mid[2]])
                window_intensity[i, j, k] = normpdf(value, 0, sigma)

    return window_intensity
    
def getSubmatrix(matrix, init, size):
    x = init[0]
    y = init[1]
    z = init[2]
    
    s_x = size[0]
    s_y = size[1]
    s_z = size[2]
    
    return matrix[x - s_x : x + s_x + 1 , y - s_y : y + s_y + 1, z - s_z : z + s_z + 1]

def conv(img, window, i, j, k):    
    window_size = window.shape
    img_size = img.shape
    w_x = window_size[0] / 2
    w_y = window_size[1] / 2
    w_z = window_size[2] / 2
	
    #Es un borde, entonces queda igual
    if (i < w_x or j < w_y or k < w_z):
        return img[i, j, k]
	
    if((img_size[0] - 1) - i < w_x or (img_size[1] - 1) - j < w_y or (img_size[2] - 1) - k < w_z):
        return img[i, j, k]
	
	#print(i, j, k)
    subimg = getSubmatrix(img, [i, j, k], [w_x, w_y, w_z])
	
    return multiply(subimg, window).sum() / window.sum()

